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ABSTRACT 

By performing local three-dimensional MHD simulations of stratified accretion disks, we investigate 
disk winds driven by MHD turbulence. Initially weak vertical magnetic fields arc effectively amplified 
by magnetorotational instability and winding due to differential rotation. Large-scale channel flows 
develop most effectively at 1.5 - 2 times the scale heights where the magnetic pressure is comparable 
to but slightly smaller than the gas pressure. The breakup of these channel flows drives structured 
disk winds by transporting the Poynting flux to the gas. These features are universally observed 
in the simulations of various initial fields. This disk wind process should play an essential role in 
the dynamical evaporation of protoplanctary disks. The breakup of channel flows also excites the 
momentum fluxes associated with Alfvenic and (magneto-)sonic waves toward the midplane, which 
possibly contribute to the sedimentation of small dust grains in protoplanctary disks. 

Subject headings: accretion, accretion disks — ISM: jets and outflows — MHD — planetary systems: 
protoplanctary disks — planetary systems: formation — turbulence 



1. INTRODUCTION 

Magnetorotation al instability (MRI; 

iBalbus fe Hawlevl 119911 ) is regarded as a robust 
mechanism to provide turbulence for an efficient out- 
ward transport of angular momentum in accretion disks. 
MHD simulations i n a local shearing box hav e been 
carried out (e.g., lHawlev. Gammie. fe Balbus I 1 19951 : 
iBrandenburg et al.|[l995l : ISano et al.ll2004|) to study the 
prope rties of MRI-driven turbulence" iMiller fc Stone I 
(|2000l ) studied vertically stratified local disks with free 
boundaries to allow leaks of mass and magnetic field. 
While their main purpose is to study general properties 
of stratified disks such as disk coronae rather than disk 
winds, they concluded that the mass flux of the outflows 
is small in the cases of initially toroidal and zero-net 
vertical flux magnetic fields. 

On the other hand, protoplanctary disks around young 
stars should have net vertical magnetic fields that are 
connected to their parental molecular clouds. In this 
case, the physical conditions of the surface of the disk 
are analogous to the open coronal holes of the Sun where 
the solar wind is driven by turbulent footpoint motions of 
the m agnetic field lines (|Sakao et al.ll2007tlTsuneta et al.l 
120081 ). Obviously, MHD turbulence excited by MRI in 
the disk is also expected to drive winds from the sur- 
faces of the accretion disk. Although such a disk wind 
mechanism may pl ay a significant role in the evolutio n 
of accretion disks (jFerreira. Dougados fc Cabrit I l2006h . 
quantitative studies have not been carried out so far be- 
cause of difficulties in the numerical treatment: a long- 
term calculation of the wind process requires accurate 
description of outgoing boundary conditions for various 
types of waves, rather than simple free boundaries 3 . 
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3 Simple free boundaries means setting the derivatives of vari- 
ables to be zero. In this case, however, unphysical reflections of 



In this Letter, we investigate disk winds driven by MRI 
with initially vertical magnetic fields utilizing the rigor- 
ous outgoing boundary condition th at was originally de- 
veloped for the simulati ons of solar (ISuzuki fc Inutsuka 1 
Eooalill) and stellar (jSuzuki 112007m winds. 

2. SETUP 

We perform 3D MHD simulations in a local shearing 
box (Hawley et al.1995), takin g into account v ertical 
stratification (jStone et al.lll996t iTurner fc Sano Il2007f) . 
We set the x-, y-, and z-coordinates as the radial, az- 
imuthal, and vertical directions, respectively. We solve 
the ideal MHD equations with an isothermal equation 
of state in a frame corotating with Kepler rotation. In 
the momentum equation we consider the vertical gravity 
by a central star, £l^z, where fio is Keplerian rotation 
frequency. We adopt a second-order Godunov-CMoCCT 
scheme, in which we solve nonlinear Ricmann problems 
with magnetic pressure at cell boundaries for compres- 
sive waves and adopt the consistent method of ch aracter- 
istics (CMoC) for the evolution of magnetic fields (|Clarkd 
[19911 . 

The simulation region is (x,y,z) = 
(±0.5H ,±2H ,±4:Ho), and is resolved by (32,64,256) 
grid points, where Hq = \/2c s /£lo is the pressure scale 
height for sound speed, c s . The shearing boundary is 
adopted for the x-direction to consider the Keplerian 
shear flow (Hawley et al. 1995). The simple periodic 
boundary is adopted for the y-direction. We prescribe 
the outflow condition in the z-dircctions by adopting 
only outgoing characteristics from all seven (six for 
isothermal gas) MHD characteristics at the z = ±4i?o 
boundaries (Suzuki & Inutsuka 2006). While the z- 
component of the magnetic flux is strictly conserved, the 
x- and y- components are not conserved because of the 
outgoing condition at the z boundaries. We initially set 

waves usually occur. For the real outgoing bounda ry, only the 
deriv atives of incoming characteristics should vanish ijThompson I 
[19871) . 
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up a hydrostatic density structure, p = po exp(— z 2 /Hq), 
with po = 1, and a constant vertical magnetic field, 
B z ,o, with the plasma /? value, (3q = 8irp c 2 /B 2 = 10 6 
(for our fiducial model) at the midplane. We use f2 = 1 
and Hq = 1, which gives c 2 s = 1/2. Small random 
perturbations, 5v = 0.005c s , are initially given for the 
seeds of MRI. 

3. RESULTS 

In most of the simulation region, the initial magnetic 
fields are moderately weak, so it is unstable with respect 
to MRI. In \z\ < 2.5Hq, however, we cannot initially re- 
solve the wavelengths, A max w 2itva/Qo-, of the most 
unstable mode because A max < Az(= Hq/32), where 
Va = B /^/AiTp is Alfvcn speed and Az is the mesh 
size (the initial A max ss 0.2Az at the midplane). First, 
MRI develops around z « ±3Hq after « 3 rotations. 
The turbulence driven by MRI gradually spreads to- 
ward the midplane, since the growth time (approximately 
cx Az/A max for Az > A max ) of the resolved wavelength 
is longer there. When t > 100 rotations, the midplane 
finally becomes turbulent. The magnetic field strength 
saturates in the entire box after t > 200 rotations, and 
the system becomes quasi-steady-state. At this time, 
Amax can be resolved even at the midplane owing to the 
increase of the field strength. We continue the simulation 
further up to 400 rotations. 

FigurcQ]is the snapshot of magnetic field (white lines), 
velocity field (arrows), and Sp/(p) (color) at t = 210 ro- 
tations, where (p) is the density averaging over each x-y 
plane and Sp = p — (p). The magnetic fields are turbu- 
lent, dominated by the toroidal (y) component because of 
winding. Angular momentum is outwardly transported 
by anisotropic stress due to the MHD turbulence. At the 
saturated state, a = {v x 5v y — B *^ p y ) /c 2 is ~ 0.01 in the 
midplane. One can also observe that the structured out- 
flows stream out from both the upper and lower bound- 
aries. Below we inspect the properties of the outflows in 
more detail. 

Figure [2] presents the mass flux of the z-component, 
pVz/poc s , in the t — z plane. We averaged pv z on the 
x-y plane at each z grid point. One can see that the 
gas flows out from both the upper and lower boundaries. 
The mass fluxes near the surface regions are highly time 
dependent with a quasi-periodic cycle of ~ 5 — 10 rota- 
tions. Moreover, from z ~ ±2Hq the mass fluxes direct 
to the midplane, almost coinciding with the periodicities 
of the outflow fluxes. In other words, the mass flows are 
ejected to both upward and downward directions from 
'injection regions' located at z ~ ±2Hq. 

These features a re consequences of the breakup of 
chann el flows (e.g., iSano fc Inutsuka 1 120011 ; iSano et all 
2004). At z ~ ±2Hq, the wavelength of the most unsta- 
ble mode with respect to MRI, A max , is comparable with 
the scale height, Hq- In the region \z\ > 2Hq, A max > Hq\ 
hence, it is stable against MRI. In the region \z\ < 2Ho, 
smaller-scale turbulence develops preferentially because 
Amax < Hq. Therefore, at z ~ ±2i?o the largest scale 
channel flows develop, and their breakup by reconnec- 
tions 4 drives the mass flows to both upward and down- 

4 We do not explicitly include the physical resistivity term in 
the calculation shown in this Letter, and so, the reconnections are 
due to the numerical effect determined by the grid scale. 



ward directions. In the region \z\ < 2Hq, the gas pressure 
largely dominates the magnetic pressure so that strong 
mass flows cannot be driven by the magnetic force asso- 
ciated with reconnections between small-scale turbulent 
fields. The periodic oscillation of 1 Kcplcrian rotation 
time around the midplane is the vertical (epicycle) mo- 
tion. 

Figure[3]presents the disk wind structure averaged over 
200 - 400 rotations. The variables are averaged on the 
x-y plane at each z point. The top panel shows that 
the average outflow velocity is nearly the sound speed 
at the upper and lower surfaces. The second panel 
presents the structures of density and plasma /3 value. 
The comparison of the final density structure (solid) with 
the initial hydrostatic structure (dotted) shows that the 
mass is loaded up to the onset regions of outflows from 
z « ±2H . In the wind region \z\ > 3H , (3 is below 
unity; the disk winds start to accelerate when the mag- 
netic pressure dominates the gas pressure. 

The third panel shows magnetic energy at the sat- 
urated state. The dashed, solid, and dotted lines 
are x-, y-, and z-components, respectively. In the y- 
component we show both mean, (B y ) 2 and fluctua- 
tion, SBy, components. {B y ) 2 is the simple average on 
the x-y planes, (B y (z)) = J j dxdyB y (x, y, z)j (L x L y ), 
and the fluctuations are determined from SB 2 (z) = 
J J dxdy(B y (x,y,z) - (B y (z))) 2 / (L x L y ), where L x (= 
Hq) and L y {= AHq) are the x and y lengths of the sim- 
ulation box. As for B x and B z the fluctuation compo- 
nents greatly dominate the means. The magnetic energy, 
which is dominated by the toroidal (y) component as a 
consequence of winding, is amplified by ks 1000 times 
of the initial value (B 2 /47r = 10~ 6 ) in most of the re- 
gion (|z| < 3Hq). While in the region near the midplane 
( | jar J < 1.5iJo), the magnetic field is dominated by the 
fluctuating component (SB y ), the mean component dom- 
inates in the regions near the surfaces (\z\ > 1.5Hq). In 
the surface regions the magnetic pressure is comparable 
to or larger than the gas pressure (J3 < 1), and so, the gas 
motions cannot control the configuration of the magnetic 
fields. Therefore, the field lines tend to be straightened 
by magnetic tension to give (B) 2 > SB 2 there, even if 
the gas is turbulent. We also note that (B 2 ) is amplified 
by MRI and Parker (1966) instability, whereas (B z ) 2 is 
strictly conserved. 

The bottom panel shows the status of energy transfer. 
The z-componcnt of the total energy flux is expressed as 

pv z Q, 2 + * + ft) + „,|L - ^(v ± B ± ), (1) 

where h is the enthalpy 5 , $ = z 2 /2, and _l_ denotes x 
and y; e.g. B\ = B 2 + B 2 . The Poynting flux is sepa- 
rated from the term of direct transport of magnetic en- 
ergy (B\v z j A.-k) and the term related to magnetic tension 
(—B z v±B±/4ir). The solid and dotted lines are respec- 
tively — B z 5vj_Bj_/4ir and -B 2 _w z /47r, where Sv± = v x and 
v y + 3/2flox (— 3/2flox is the background Kepler rota- 
tion). The dashed line is the potential energy term, pv z &. 
The dot-dashed line is the energy flux of sound waves, 
Sp8v z c 2 (see below). The gas pressure (pv z ) and hy- 
drodynamical turbulent pressure (pSv 2 v z /2) are smaller 

5 Formally, h = i? s log p for isothermal gas 




Fig. 1. — Snapshot of the local disk structure at t = 210 rotations. The white solid lines are magnetic fields, the arrows indicate velocity 
fields, and the colored region corresponds to Sp/ (p) > 0. 
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Fig. 2. — Time-height diagram of the mass flux, pv z , normalized 
by pqc s - pv z is averaged on the x-y plane at each height, z (vertical 
axis). The unit of of horizontal axis is the rotation period (2tt/Qq). 



than these terms. The kinetic energy flux of the winds 
(^pv z ) is also small < 10 -5 at the outer boundaries. 

The figure shows that the materials near the surfaces 
are lifted up by the conversion of the Poynting flux; the 
absolute values of the Poynting flux terms (solid and dot- 
ted) decrease with height in the region \z\ > 2Hq, and the 
absolute value of the potential energy flux (dashed) in- 
creases. Both magnetic pressure and tension terms con- 
tribute almost equally. 

—B z 6v±B±/4tt and 8p8v z c 2 are the net energy fluxes 
of Alfven waves and sound waves to the +z-dircction. 
—B z Sv±B±/4:TT can be rewritten as 



— B z Sv ± B± = pv A , z (Sv ± , - 5v ± _), 

47T 



(2) 



where da. 2 = B z /^/4irp, &nd8v± t ± = ^ (8v± ^fB± / ^/Anp) 
are Elsasser variables, which correspond to the ampli- 
tudes of Alfven waves propagating to the ±z-directions. 
Sp8v z c 2 is also rewritten as 



where 5v\\ 



8p8v z c 2 s = pc s (8v 2 + - 8v 2 _), (3) 

; ||.± = \{& v z ± °s^f) denote the amplitudes of 
sound waves 6 propagating in the ±z-directions. 

An interesting feature of —B z 8v±B±/4ir is that the 
sign changes at z w ±1.6i?o (the circles in the bot- 
tom panel); in z > 1.6Hq(< — 1.6Ho), the flux is up- 
ward (downward) and \z\ < 1.6Hq the flux is toward 

6 Strictly speaking, these are magnetosonic waves, namely the 
fast mode in the high /3 plasma, and the slow mode that propagates 
along 2 in the low /3 plasma. Note also that the signs are oppo- 
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Fig. 3. — Time-averaged disk structure during i = 200-400 ro- 
tations. The variables are also averaged on the x-y plane at each 
z grid. The top panel shows v z /c a (solid), whereas the dotted line 
is the initial condition (v z /c s = 0). The second panel presents 
density (solid; left axis) and plasma f3 (dashed; right axis), in 
comparison with the initial condition (dotted; for both density 
and plasma /3). The third panel presents the magnetic energy, 
B 2 /Atv. The dashed, solid, and dotted lines correspond to the 
x-, y-, and ^-components, and the y-component shows both mean 
(thick) and fluctuation (thin) components. The bottom panel illus- 
trates the energy flux in units of po(-Ho^) 3 ■ The solid and dotted 
lines are the Poynting flux associated with the magnetic tension 
(—B z 5vj_B^/4tt) and the magnetic energy (B^v z /4-7r). The dot- 
dashed line is the net energy flux due to sound waves (5p&v z c 2 s : 
see the text). The dashed line is the term concerning the potential 
energy (pv z Q). The circles are the 'injection regions' defined as the 
locations where the signs of — B z 5v±B±/4n change. 
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the midplane. This is a consequence of the breakup of 
channel flows, as described previously. Reconnections 
break up channel flows and generate large amplitude 
(|<5-Bj_| > \B Z \) Alfven(ic) waves in both upward and 
downward directions. Sp8v z c 2 . is also directed to the mid- 
plane, namely sound(-like) waves are generated by the re- 
connections. The absolute value of the energy flux peaks 
at the injection regions at z « ±1.6ffo- The energy flux 
of sound waves to the midplane is larger than that of 
the Alfven waves. On the other hand, the Alfvcn wave 
component greatly dominates in the flux to the surfaces. 
This is because (3 > 1 (the gas pressure dominates) in 
\z\ < 2Hq and j3 < 1 near the surfaces. 




Fig. 4. — Dependence of the final (3 structure on initial magnetic 
fields. The dash-dot-dotted, dot-dashed, dotted, solid, and dashed 
lines correspond to the initial vertical fields with /3o = 10 4 , 10 5 , 10 6 , 
and 10 7 , and the initially toroidal field cases. The arrows indicate 
the locations of the injection regions defined as the points where the 
sings of —B z 5v±B± change. Longer arrows correspond to smaller 
Po (the shortest arrows are for the initially toroidal case). 



In order to study the effects of the initial magnetic 
fields, we performed the simulations with initial different 
vertical field strengths, f3 = 10 4 , 10 5 , 10 6 , 10 7 , and ini- 
tially toroidal field in \z\ < 3Hq with (3q = 10 6 at the 
midplanes. Figure [4] compares {(3) structures. Figured] 
compares the sum of disk mass fluxes from the upper 
and lower boundaries with various cases for the initial (3q 
at the midplane. In all cases except the Pq = 10 4 case, 
we averaged the variables during 200 rotations after the 
quasi-steady-states are achieved. For the /3o = 10 4 case, 
we show the time averages during 25-55 rotations be- 
cause the quasi-steady-state is not achieved but 90 % of 
the total mass escapes at 75 rotations as a result of the 
effective mass loss by disk winds. 

All the cases of (3q > 10 5 show very similar ((3) struc- 
tures (Figure [4|): (3 w 100 at the midplanes, namely the 
magnetic energy can be amplified to « 1% of the gas 
energy. The f3 values decrease with increasing height 
mainly because of the decrease of the densities (the mag- 
netic energies stay nearly constant; Figure [3]). When 
the magnetic energy dominates (/3 < 1), the disk winds 



o r-o 

\ o 



3 * 

^ N o 



o- 



O 



O 



: a 



o 



10' 



Initially 
Toroida 



10° 

Initial 8 
Constant B_ 



10- 



10" 



Fig. 5. — Sum of the mass fluxes normalized by poc s of the disk 
winds from the top and bottom boundaries of the simulation box. 
The horizontal axis indicates the initial /3o at the midplane for the 
vertical field cases circles. The initially toroidal case is plotted at 
the leftmost location (triangle). 



start to accelerate. In the wind regions, the (3 values 
stay at (3 = 0.1-1, rather than further decrease, owing 
to the increased density by the lift-up gas in the winds. 
The locations of the injection regions also concentrate at 
\z\/Hq = 1.5-2 except for the initial (3q = 10 4 case. The 
plasma (3 values of the injection regions are 1-10; the 
magnetic energy is comparable to but slightly smaller 
than the equipartition value. This condition is favorable 
for driving mass motions by the breakup of large-scale 
channel flows; if the magnetic energy is larger, the field 
configuration becomes more coherent due to the tension 
so that reconnections hardly occur; if the magnetic en- 
ergy is smaller, the reconnections cannot drive strong 
mass motion. The mass flux of the disk winds (Figure 
[5]) only weakly depends on (3q for f3o > 10 6 , while it in- 
creases for smaller (3q almost in proportion with /3q. 

4. DISCUSSIONS 

We have shown that the gas is lifted up from the in- 
jection regions at \z\/H w 1.5-2 to the surfaces by the 
Poynting flux and streams out from the upper and lower 
boundaries of the simulation box. In effect, our results 
determine the condition of "mass loading" in various 
models of global disk wind ( e.g., iBlandford fc Pavnel 
fl98l iKudoh fc Shibata I H9971 Ferreira et al 2006), in 
which one usually fixes the mass flux in advance by set- 
ting the densities at the 'bases' of winds. With our outgo- 
ing boundary condition, we implicitly assume that once 
the gas goes out of the z-boundaries of the simulation box 
it does not return. The validity of this treatment needs 
to be examined by the global modeling of accretion disks. 

Although in the shearing box treatment the time- 
averaged net mass flow in the x-direction is zero, we can 
estimate the accretion velocity, — v r ~ ac 2 s /rQ,Q, from 
the angular momentum balance under steady states, 
where r is a cylindrical distance from a central star 
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(|Shakura fe Sunvaev 1119731 ). The ratio of the mass-loss 
rate, M z , from the simulation box by the disk winds to 
the mass accretion rate, M r , passing through the y-z 
plane becomes 



(Ms) J J dxdy(pv z ) _ {{pv z ))L x rVl a 



(M r ) Jdy(m)c 2 s /(rfl ) 
lOHoJ \H 



«Ea»c 



:0.05 



(4) 



where E = J dzp and a = J dzpa/T,; () denotes the time- 
average and (()) is the average of time and x-y planes. 
The final value is estimated from ((pv z )) w 8 x 10 _5 poC s 
and a ~ 0.012 in our fiducial run for a moderately thin 
disk with Hq/t = 0.1. We should take the above estimate 
as an upper limit because the calculation with a larger 
vertical box size may give a smaller mass flux at the 
top and bottom boundaries. The angular momentum 
loss rates in vertical and radial directions give the same 
scaling, 

(C z ) _ J J dxdy(pv z )r 2 n Q _ ((pv z ))L x rn (M z ) 



(Cr 



Jdyr(T,a)c 2 s 



«Ea»c 



(Mr) ' 
(5) 

because the time-averaged specific angular momentum 
carried by the winds is approximately the same as that 
in the disk material at the same radial position in the 
shearing box treatment. In a realistic situation, however, 
winds possibly carry a larger specific angular momentum 
than the disk material. For such studies, we need to 
model global accretion disks with disk winds, which is 
also important from the viewpoint of an gular momentum 
evolu tion of the star-disk system (e.g. iMatt fc Pudritz 1 
l2005h . 

Hereafter we discuss the evolution of protoplanetary 
disks, as an application of our results. As a reference 
model, we use the minimum-mass solar nebula (MMSN) 
of Hayashi (1981), which gives the midplane density, po = 

1.4 x 10 -9 (jxa) U g cm~ 3 . Then, the initial vertical 
magnetic field of (3q — 10 6 in our fiducial run corresponds 
to B z o ~ 0.01 G, and the saturated field strength is 
B w 1 G at 1 AU. 

First, we examine how much the disk wind contributes 
to the evaporation of protoplanetary disks (see e.g. 



iDullemond et al.1l2007l for other mechanisms). After the 
saturation of the magnetic fields, f» 5% of the total disk 
mass is lost from the simulation box from 200 to 400 ro- 
tations by the disk winds in our fiducial case. Assuming 
a disk around a central star with the solar mass (1 ro- 
tation = 1 yr at 1 AU), we have the timescales of the 
evaporation, t cv w 4000 yr at 1 AU, and 6 x 10 5 yr at 
30 AU. Although this is rather short in comparison with 
recent observational results (typic ally r ov ~ 10 6 " 7 yrs, 
e.g., lHaisch. Lada. fc Lada 1120011) . this is not a severe 
contradiction because we have not yet taken into account 
the global radial accretion of the disk mass, which contin- 
uously supplies the mass from the outer region. Another 
important issue that affects, and might reduce, the mass 
flux of disk winds is the effect of resistivity, which re- 
quires an additional detailed analysis of the ioniz ation 
structure (|Sano et al.ll200Ct llnutsuka fc Sano Il2005t ) and 
will be the scope of our next paper. Here, the estimated 
r ev should be taken as a lower limit. 

The disk scale height has a relation of Hq/t oc r 1 / 4 for 
the MMSN. Combining with Equation (|2|), we infer that 
the dynamical evaporation by disk winds, in comparison 
with accretion, becomes relatively more important in the 
inner parts of protoplanetary disks than in the outer re- 
gions for a constant initial /3o structure (B z0 oc r~ n / 4 
for the MMSN). 

Finally, we should point out the effects of waves on 
dusts in protoplanetary disks. We have shown that the 
momentum flux of Alfvcnic and sound-like waves directs 
to the midplane from the injection regions. The momen- 
tum flux of the sound- like waves (Spdv z ) can push dust 
grains to the midplane by gas-dust collisions. Dusts are 
usually weakly charged; in this case Alfvenic waves also 
contribute to the sedimentation of dusts to the midplane 
through ponderomotive force or du st-cyclotron resonance 
(|Vidotto fc Jatenco-Pereirall20060 . 
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